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Abstract 

We introduce an implicit solvent Molecular Dynamics approach for calculating ionic fluxes 
through narrow nano-pores and transmembrane channels. The method relies on a dual-control- 
volume grand-canonical molecular dynamics (DCV-GCMD) simulation and the analytical solution 
for the electrostatic potential inside a cylindrical nano-pore recently obtained by Levin [Europhys. 
Lett., 76, 163 (2006)]. The theory is used to calculate the ionic fluxes through an artificial trans- 
membrane channel which mimics the antibacterial gramicidin A channel. Both current- volt age 
and current-concentration relations are calculated under various experimental conditions. We 
show that our results are comparable to the characteristics associated to the gramicidin A pore, 
specially the existence of two binding sites inside the pore and the observed saturation in the 
current-concentration profiles. 

PACS numbers: 87.16.Uv, 87.10.Tf, 87.16.A- 
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I. INTRODUCTION 



Ion channels are structures formed when specific proteins are incorporated into the phos- 
pholipid membrane The channels serve to establish an electrostatic potential gradient 
across the cell membrane by allowing an ion specific flux to pass through the membrane. 
There are many different ion channels in living cells. They differ in composition, pore struc- 
ture, and ion selectivity 

a a. 

Thus, a full description of the architecture and operation 
of ion channels is a very difficult task. From purely electrostatic point of view, operation 
of ion channels presents an interesting theoretical puzzle. Since the channel passes through 
a low-dielectric membrane, there exists a large potential energy barrier for ionic solvation 
inside the pore. Yet, in practice, it is well known that when open, ion channels sustain a 
very large ionic transport rate, compatible with a free diffusion. 

A recent advance in nanoscale technology is the synthesizing of nanotubes with chemical 
modified surface [41, which allows to construct functionalized nanotubes that mimic the 

fin 

behavior of biological ion channels and are stables when inserted in a lipid bilayer |5|, |6j, and 
also exhibits a high energy barrier [7]. An excellent review of recent theoretical advances in 
synthetic nanotubes which mimic the behavior of biological ion channels is given by Hilder et 
al fl. 

There have been a number of different strategies proposed to understand the ion trans- 
port across biological and synthetic channels. One approach used extensively over the last 
few years is all-atom molecular dynamics simulation (MD). The advantage of atomistic MD 
is that molecular structure of the pore, ionic species, and water are taken explicitly into ac- 
count {2I, [7 -[l4]. Although this method is very appealing, it remains a huge task to relate the 
measured observables to the experimental results. One of the main obstacles are the compu- 
tational times needed to achieve the time scales observed experimentally |15l- 18] which, in 

19l | . Another 



most cases, requires large-scale MD simulations in special-purpose machines 
problem is the choice of the molecular force field to be used. For instance, in the case of the 
biological channel gramicidin A (gA), one of the most studied ionic channels, the first MD 
simulations have predicted a potential of mean-force (PMF) with a barrier much larger than 
expected experimentally 1^]. Recent MD simulations with molecular force field properly 



constructed appear to produce PMFs in semiquantitative agreement with experiments 
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22] . Although promising, this agreement seems to be very sensitive to the molecular force 
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field used and with the membrane where the gA channel is embedded 22J. Furthermore, 
the classical water models used in all-atom simulations are parametrized for bulk water 
and might show erroneous behavior in a strongly confined environment. Such artifacts of 
classical water and ion models have been recently observed in the studies of ionic solvation 
near interfaces, when compared to the full ab initio simulations 23). It has been shown 
that for ionic solvation in an interfacial geometry properly constructed dielectric continuum 
models agree better with results of full ab initio simulations than the classical explicit water 



models 



23 



One alternative to atomistic MDs are the, so called, Brownian dynamics simulation 



(BD) 10|, |26|, [27| . In these simulations only ionic movement is integrated, while the protein 



degrees of freedom are held fixed and the water is treated as a uniform dielectric continuum. 
This significantly reduces the computational cost of the simulation, allowing to access much 
larger time scale. Nevertheless, BD simulations still requires solution of Poisson partial- 
differential equation at each new time step of the simulation, making them quite difficult 
to implement. Recently new atomic-resolution BD simulations have been proposed 28]. 
Using PMFs derived from all-atom MD simulations, this method was able to access large 



28]. 



simulation times at reasonably low computational cost 

Thus, in order to avoid the difficulties of all-atom and BD simulations, as well as to be able 
to explore large time scales necessary for measuring the transmembrane currents, it might be 
useful to use the dielectric continuum models of ion channels as a first order approximation 
for the transmembrane dynamics. In this paper we explore the transmembrane ion fluxes 
using the continuum electrostatics model recently introduced by Levin [29]. Briefly stated, 
Levin solved the Poisson equation with the appropriate boundary conditions to obtain an 
analytical expression for the electrostatic interaction potential (Green function) between the 
charges inside a finite cylindrical pore passing through a low-dielectric phosphoric membrane. 
This electrostatic potential can be used to calculate the forces acting between the ions 
inside the pore and between the ions and the charged protein residues embedded in the cell 
membrane. It is our goal in this study to use the interaction potential derived by Levin in a 



series of Dual-Control- Volume Grand- Canonical Molecular Dynamics (DCV-GCMD) 30l.l31| 



simulations of a simple model of narrow transmembrane channel. We should stress t 
our model is very different from the mean-field Poisson-Nernst-Plank (PNP) theory 



hat 
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33] . where the Poisson equation and the continuity equations for mobile ions are solved 
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FIG. 1: Schematic depiction of the simulation cell. The ionic channel is built as a cylinder made of 
a sequence of annular rings of LJ spheres. The system is confined in x direction by two flat walls 
(represented by gray) at each end of the simulation box. The same flat walls are used to maintain 
the flux of cations between the two reservoirs (CVs) through the channel. Periodic boundary 
conditions are applied in both y and z directions. 

simultaneously in a self-consistent way. For narrow channel such as gA the correlation 
effects between the ions are of fundamental importance and the mean-field description of 



ionic conduction is bound to fail 



34] 



Our simulations are performed for different electrolyte concentrations and external applied 
electric fields to verify if our model reproduces the ionic flow obtained for one specific system: 
the gramicidin A channel. The comparision between our model and experiments for various 
concentrations and external fields is an indication that it is able to reproduce the observed 
ionic fluxes in narrow channels Q, . 

The paper is organized as follows. The model and computational details are given in 
Sec. [TT1 The results are discussed in Sec. [Till and the summary and conclusions are presented 
in Sec. HV1 

II. THE MODEL SYSTEM AND THE SIMULATION METHODOLOGY 
A. A model for transmembrane channels 

We use the channel-reservoir setup illustrated schematically in Fig. [TJto calculate the ionic 
currents due to electrostatic potential gradients. The simulation box, a cubic parallelepiped 
with volume 5LxLxL and L = 20 A, contains two reservoirs (two control volumes CVi and 
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CV2, within which the chemical potentials of the ionic species are maintained constant) and 
a membrane with a channel connecting these two reservoirs through small buffers regions. 
The channel structure, a simplified model of a transmembrane pore, is built as a cylindrical 
tube, with radius a = 3 A and length L c = 35 A, made of stationary Lennard- Jones (LJ) 
spheres of diameter o c = 2 A. 

Both sides of the channel structure — except for the orifices — are bound by confining 
walls, as well as the two extremes of the simulation box in the x direction, see Fig. [TJ 
The system contains positive and negative ions with diameters cr+ and <r_, respectively, 
inside a structureless solvent of dielectric constant e w = 80 — the same value of dielectric 
is used inside and outside the channel — while the membrane has a dielectric constant 
equals to e m = 2, both in units of vacuum permittivity. It should be noted that in confined 
environments the dielectric constant of water could be considerably lower than the bulk value 
80 0, [37)]. However, in this work we follow the same prescription of previous dielectric 
continuum models, where the use of same dielectric constant is compensated by fixing a 



lower diffusion constant for the positive ions moving inside the channel [10 . 

The particle-particle interactions are separated into short and long-range contributions, 
while the particle-channel has only a short-range interaction. For the short-range part we 
will use the WCA LJ potential |39| 



U L j{r) - U LJ (r c ) , r<r c , 
, r > r. 



U% CA (r)= { LJVJ LJVCJ ' " c ' (1) 



C ) 



where U^j(r) is the standard 12-6 LJ potential. The cut-off distance is r c = 2 1/f V i7 -, where 
&ij = (&i + 0j)/2 is the center-to-center separation between an ion of species i (cation 
or anion) and a particle of species j (cation, anion, or a fixed LJ channel sphere) sepa- 
rated by a distance r. The confining walls in simulation box extremes and surrounding 
the channel structure are modeled with the same WCA LJ potential, however considering 
the x-projection of the distance between one ion in the bulk and the wall position. The 
long-range contribution is calculated depending on the region where the ion is located. For 
the regions outside the channel, the interaction energy between the two ions is the usual 
Coulomb potential 



= - (2) 

where r« is the distance between the two ions. The infinite extent of the particle reservoir 



is taken into account using the Ewald summation. 

For the region inside the channel we will use the model introduced by Levin [29]. In 
this model ions inside the pore interact through the electrostatic potential which is only a 
function of their separation in the x-direction (along the axis of symmetry of the channel). 
This is quite reasonable for narrow channels that allows only single file motion. When one 
ion with charge q enters into the channel it interacts with the other ions and with the p 
residues of charge — q embedded into the channel wall at transverse distance p from the 
central axis. In addition to the electrostatic interaction between the ions and the residues, 
there is a self-energy penalty U associated with an ion entering into the region surrounded 
by the low-dielectric material. The electrostatic energy of the ions inside the pore is then 

^ N N N v N 

V = fcVtofai, x i) +^2^2^0X1^, p,xj) + S ^q l U(x i ) , (3) 

i=l j=l i=l 3=1 i=l 

where (fi n (xi, Xj) = (pifa, Xj) + ^{xi, Xj). The two electrostatic potentials are {22 1 

<fl(Xi, Xj) = 

q f°° {a 2 (k)e k ^- x ^~ 2kLc + a{k)P{k)[e~ k ( Xi+x ^ + e k ^+ x ^- 2kL ^] + (3 2 (k)e~ k ^- x ^} 

(3 2 (k) -a 2 (k)exp[-2kL c ] ' 

(4) 
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where a = [k — [k 2 + K 2 ) 1 / 2 ] and (3 = [k + [k 2 + k 2 ) 1 ' 2 ], with k the inverse Debye length that 
characterizes the electrolyte concentration in the two reservoirs, and 

, x s _ 4g(e w - £p) KotknajK^knCi) sm(k n Xj) sm(k n Xj) 
2 " 1 e w L c £^e w Ii(k n a)K (k n a) +e p I (k n a)Ki(k n a) 7 1 

where J n and K n are the modified Bessel functions of order n. Here a is the channel radius 
and k n = tvk/L c . The interaction potential between an ion and a fixed charged residue is 

, p x x 4? y> JL K (k n p) sm(k n Xj) sm(k n Xj) 

Finally, the electrostatic potential responsible for the barrier is given by 

_q_ r {2a 2 (k)e- 2kL c + a(k)(3(k)[ e - 2kx + e 2k{x - L ^}} 
{X) ~2^J P 2 (k)-a 2 (k)exp[-2kL c ] + 



2q(e w - ep) y> jjf (fc n o)A"i(fc n a) sin 2 (£; ra x) ^ qn_ 

, ewh(k n a)K (k n a) + t p Io(k n a)Ki(k n a) 2e w ' 



(k n a)K(j(k n a) + e p Io(k n a)Ki(k n a) 2e v 
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In addition we will place two residues of charge — q each embedded into the channel wall 
(p = 3 A) at positions x = —10.5 A and x = 10.5 A, to represent the two binding sites of 
our transmembrane channel. The residues are placed inside the channel wall leaving space 
for the mobile ions to flow. The channel length and radius are comparable with the size of 
nanotubes and other models for the biological channels such as gramicidin A (gA) Q, 



Although the 



groups 17|, |40j, these behave similarly to the charged binding sites used in our model [lOj. 
Furthermore, recent studies have shown that it is possible to use functionalized synthetic 
nanotubes to obtain ionic fluxes similar to the ones observed in gA channel |35 ] - 

B. Simulation details 

Ion channels can be highly selective to the cation passage, and nanotubes can be func- 
tionalized to mimic this behavior l|, |8j. Although the actual transport process is very 
complicated, we can identify three main steps (1) cation entry, where the positive ions 
are dehydrated, (2) cation translocation through the channel, and (3) cation exit. Since in 
our model we do not consider water molecules explicitly, we simulate step (1) above using 
a diameter for the cations equals to <j + = 2 A, while the negative ions have a diameter 
cr_ = 4 A. Since the channel is very narrow, the available radius for ion movement is ap- 
proximately a — a c /2 in Fig. [TJ Therefore, using these diameters only positive ions can enter 
into the channel. In all interactions the energy parameter of LJ potential was defined as 
e = lksT. During the MD steps, we have used Langevin dynamics to simulate the effect of 
solvent on the cation and anion movement, solving the equation of motion of ion i, 

rrii Si = Fi — m^Vi + W(t) , (8) 

where Fi is the total force on ion % due to all entities explicitly present in the model (other 
ions, protein residues, and walls), 7 is the friction coefficient, and W(t) is the random 



force 



41] due to solvent. The temperature of the system is maintained constant using the 



fluctuation-dissipation theorem, 

(W(t).W(t')) = 6k B T^5(t - t') , (9) 

which relates the friction coefficient to the fluctuations of the random force using the appro- 
priate 7 value. In our simulations we have used m,i = 6.5 x 10" 26 kg, corresponding to the 



mass of K + ion; 7 _1 = 53 fs for the region outside the channel, corresponding to diffusion 
constant D = kBTjmcj = 3.37 x 10~ 9 m 2 /s; and 7" 1 = 7 _1 /3 inside the channel corre- 
sponding to diffusion constant of D — 1.12 x 10 _9 m 2 /s. These values were chosen in order 
to maintain the temperature fixed at 298 K and to reproduce the experimental behavior, 
particularly the saturation observed in the current-concentration curves [10| . 

We apply a linear voltage gradient across the membrane from the right border of CVi to 
the left border of CV2. The concentrations in both reservoirs are maintained constant using 



DCV-GCMD simulations 



30 



31 



42 



4J|. Briefly stated, in the DCV-GCMD simulation 
two control volumes (CVs) are initially prepared at desired concentrations using the grand 
canonical Monte Carlo (GCMC) simulation and then evolved in time using the molecular 
dynamics (MD) simulation. Since the dynamics alters the CV concentrations, the MD steps 
are intercalated with the grand canonical Monte Carlo (GCMC) simulations performed inside 
the two control volumes (CV) shown in Fig.[TJ This restores the concentrations to their initial 
values. In our simulations we have used 50 GCMC steps for every 500 MD steps. It should be 
noted that our DCV-GCMD method is closely related to the GCMC/BD algorithm proposed 
by Im et al. [47] to simulate the ionic conductance in membrane channels. 

Since our simulation setup is periodic only in y and z directions, as shown in Fig. Q], for 
the long-range interactions described by Eq. (j2j) we have used the Ewald summation with the 



implementation of Yeh and Berkowitz 48| for the slab geometry. The equations of motion 
were integrated using the velocity Verlet algorithm, with a time step of 8.0 fs in the MD 
simulations. The observables were obtained using 5 to 10 independent realizations. 

III. RESULTS AND DISCUSSION 

We start our discussion with the two CVs having the same concentration, and the ionic dif- 
fusion through the channel driven by the externally imposed electrostatic potential gradient 
between the two CVs. We are particularly interested in the current-voltage and current- 
concentration curves and their comparison with the experimental results, mainly the ap- 
pearance of the experimentally observed saturation in the current-concentration profiles. In 
Fig. [2] we show the current-voltage curve for the concentration of 0.5 M in both CVs. As 
one can see, we obtain the same expected linear dependence between the ionic current and 



applied voltage, observed in experiments and earlier BD simulations for gA channel lOj 
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FIG. 2: Current-voltage curve for 0.5 M concentration of electrolyte in both CVs. Open circles 
are our simulations results . The open squares and the open triangle are the experimental and BD 
results for gA channel, respectively, extracted from Fig. 12 of Ref. [lo| . 

Next we analyze the behavior of the current-concentration profiles for two externally 
100 mV and 200 mV, Fig. [3] Compared to the experimental and BD 



applied voltages o 

simulation results 10j, our model is able to capture the same saturation of the ionic current. 
For 200 mV and above 0.9 M concentration, however, we find some deviation from the 
experimental results. The deviation appears to be an artifact of including only 30 terms in 
the infinite series of the electrostatic potential inside the cell, Eqs. ([5]) and (J6]). Cutting 
the infinite series at only 30 terms softens the repulsion between the ions inside the channel 
favoring an increased flux. 



Experimental [49M52| and BD simulations [10| data for gA channel and theoretical re- 
sults for functionalized nanotubes [ssj] have proposed two large concentration peaks at the 
binding sites separated by a cation depleted region. This is exactly what we observe in our 
simulations, as shown in Fig. H]for 0.5 M monovalent solution in both CVs with no voltage 
and 200 mV applied potential difference. 

To understand better the mechanism of ionic translocation, in Fig. [5] we plot the elec- 
trostatic potential inside the transmembrane channel, Eq. ([3]), at zero applied voltage and 
for different occupancy inside the channel and residues embedded in the channel structure. 
These profiles were obtained using a dynamical PMF, where an ion, called ion*, is pulled 
through the channel, moving a displacement Sx each time. If there are other ions inside the 
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FIG. 3: Current-concentration curves for external voltages of 100 mV and 200 mV. Our simulations 



are represented by open circles. The open squares and 



results, respectively, extracted from Fig. 12 of Ref. 



open triangles are the experimental and BD 
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FIG. 4: Mean number of cations in the axial direction of gA channel with (a) no applied voltage 
and (b) 200 mV applied potential. In both cases there is 0.5 M monovalent electrolyte solution in 
the two CVs. 

channel (called ions 1, 2, 3, etc) we move the ion*, and then we allow the system to relax 
for 4 ps. During this relaxation time only ions 1, 2, etc. can move, while the ion* remains 
fixed. Then we evaluate the mean force on ion*, and perform the next displacement 5x. 
The procedure is repeated until ion* leaves the channel. Note that during this process the 
other ions inside the channel will feel the electrostatic repulsion from ion* (recall that only 
cations are allowed inside the channel), and will also be forced out of the channel. 
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FIG. 5: Electrostatic energy profile inside the channel obtained from Levin's model [29]. The curves 
show the potential of mean force (PMF) felt by an ion moving through a gA-like channel with two 
charged residues at x = —10.5 A and x = +10.5 A (solid line). The PMF inside a channel which 
already contains a free cation near the first binding site is represented by a dashed line, while the 
PMF when there are two free cations near the two binding sites is represented by a doted line. The 
dot-dashed line represents the electrostatic potential produced by Eq. ([7J). The concentrations in 
the two CVs are 0.5 M, with no applied voltage. 

If the channel has no fixed residues and is empty, equation ([7j) leads to a large electrostatic 
potential energy barrier of approximately 8 ksT, which prevents any cation entrance. On 
the other hand, if the protein has charged residues embedded into the surface of the channel, 
the scenario changes completely. As one can see in solid line in Fig. [U with no other cation 
inside the channel, the ion* feels an energy barrier of approximately 4/c#T separating two 
deep wells at the positions of the binding sites. 

Entrance of a cation alters drastically the potential energy landscape. Suppose that one 
ion, ion 1, is already inside the channel. We are interested in the potential of mean- force 
(PMF) that a second ion, ion*, will feel as it moves through the channel. This PMF is 
plotted as a dashed curve in Fig. [5j A short time after the ion 1 has entered the channel 
its most probable location is at the first binding site. Thus, when ion* enters the channel, 
it sees the field of attractive residue partially screened by the ion 1, so that the depth of 
the potential well produced by the first residue is significantly smaller. As the ion* moves 
farther into the channel, it forces the ion 1 to move to the position of the second residue 
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and, eventually, to completely leave the channel. Consider now that there are two ions, 
ion 1 and ion 2, inside the channel, with ion 1 located at the first binding site and ion 2 
at the second binding site. The ion* will then encounter a flat energy landscape shown in 
Fig. [5] by the dotted curve. This then explains the fast transport of ions through the ion 
channel observed experimentally — the cations present at the binding sites screen both the 
electric field produced by the residues and by the induced charge on the channel wall. In 
the BD studies [lo| the saturation was observed as the result of an imposed double well 
potential of depth of 8 ksT and a barrier between the wells of 5 ksT. In the present study, 
the electrostatic energy landscape was obtained self-consistently, assuming the existence of 
two monovalen binding sites inside the channel. 

In Fig. [6] we show the distribution of times of cation translocation through the channel 
for voltage difference of 100 mV and CV concentration of 0.5 M. The mean first passage 
time (MFPT) as a function of applied voltage and CVs concentration is also shown in Fig. |6j 
The saturation observed in Fig. |3] can also be seen in the MFPT results. 
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FIG. 6: (a) Distribution of passage time Pit) for 100 mV applied voltage and CV concentration of 
0.5 M. (b) Mean first passage time (MFPT) as a function of CV concentrations for 100 mV (open 
circles) and 200 mV (open squares) applied voltages. 



IV. CONCLUSIONS 

We have used a Dual-Control- Volume Grand- Canonical Molecular Dynamics (DCV- 
GCMD) simulations to study the flow of ions through narrow pores across low dielectric 
membranes. To account for the electrostatic interactions inside the channel we have used 
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the analytical potential recently derived by Levin |29j. For electrolyte concentrations up 
to 1 M our model is able to reproduce most of experimental behavior of the ionic current 
as a function of electrostatic potential difference. To obtain reliable results for larger con- 
centrations one must include more terms in the infinite series for the ion-ion interaction 
potential. For physiological concentrations of electrolyte, on the other hand, the continuum 
model presented in this work appears to show excellent results. The utility of the model is 
that the simulations do not require expensive computational resources and can be run on a 
desktop PC. One interesting application of the model is to study the dependence of ionic 
current on mutations of charged residues — the position and the charge of the residues that 
enter into the electrostatic potential can be optimized to obtain the desirable channel char- 
acteristics. This could be of interest for design and implementation of molecular engineered 
of ion channels and nano-pores, biological or synthetic. 
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